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We present a new method of ab-initio electronic- 
structure calculation including the spin-fluctuation (SF) self- 
consistently. We start from the Luttinger-Ward functional 
given as the sum of the LDA functional plus the temperature- 
dependent part of the SF energy functional. The size of inter- 
actions used in it are determined in a similar manner to the 
self-consistent renormalization theory by Moriya and Kawa- 
bata. Obtained paramagnetic susceptibilities on Pd, Ni, Fe, 
and fcc-Co above T c show rather good agreements with ex- 
periments. 

In magnetic metals like Fe, LDA gives quite good de- 
scription for ground state properties at zero-temperature 
T = 0; LDA reproduces the so-called Slater-Pauling 
curve obtained by experiments very well [1], also does 
the spin waves [2] . But neither LDA nor its crude exten- 
sion to finite temperature contains the contributions due 
to the spin- fluctuation (SF). SF for magnetic materials 
should have strong temperature dependence, and enforce 
them paramagnetic above their critical temperatures T c . 
In Refs. [3] and [4], methods including SF had developed 
with adiabatic approximation in an ab-initio scheme; 
they first calculated the adiabatic surface of all the con- 
figuration of static SF, then calculated the partition func- 
tion assuming spin variables as classical. It works rather 
well, though these approximations could be problematic 
in paramagnetic phases at high T where the Stonar-like 
excitations and spin-wave like excitation could mix well. 
In addition, the ground state, where they calculated the 
adiabatic surface, can be changed at higher T. Especially 
in Pd, which is "almost-ferromagnetic", these methods 
have not been applied to, presumably because the change 
of ground state would cause problems. As the oppo- 
site limit where the Heisenberg-like model works well, 
we know "metallic weak- ferro- magnet" , which are well 
described by the self-consistent renormalization (SCR) 
theory [5] given by Moriya and Kawabata. They evalu- 
ate SF in a kind of random-phase- approximation (RPA), 
where we should determine the size of interaction in a 
self-consistent (SC) manner. Here we present a new ab- 
initio SC method including SF, where the size of interac- 
tion is determined self-consistently as in the SCR theory. 
As results compared with experiments, we show param- 
agnetic susceptibilities for Pd, Ni, Fe, and fcc-Co above 
T c . Our calculations are based on a finite-temperature 
version of AkaiKKR [6], in the atomic sphere approxi- 
mation (ASA). 



In the treatment like LDA+[7 and their extensions 
like LDA+ 'dynamical mean field theory (DMFT)', a key 
problem is in the determination of the effective inter- 
action U. Lichtcnstcin, Katsnelson, and Kotliar pre- 
sented such a scheme to describe magnetic susceptibil- 
ities of Ni and Fe at Hight T. Their method can take 
into account higher-level of correlation effects than our 
treatment, though their U is given externally [8]. Very 
recently, Bicrmann, Aryasetiawan, and Georges sug- 
gested GVF+DMFT to determine them self-consistently 
[7] . Our method can be identified as a simplified version 
of such a method. 

We start from the Luttinger-Ward (LW) functional 
£l[g] at finite temperature. It can be written as 
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Here g a (l,2) is the finite-temperature Green function, 
and 1 = riTi. n a (r) denotes the electron density and 
is treated as the functional of g a (l,2) in Eq.(l). We 
may identify the first term in right-had side (RHS) in 
Eq.(l) as "the population-entropy" term because it is 
reduced to be — T times the entropy of the system in 
the non-interacting case. The second term is the kinetic 
energy term, and i? ext [n] is the external-potential term 
containing the chemical potential. Terms Eqou and E KC 
are related to the Coulomb interaction between electrons. 
Sum of them are referred to as $ [g] in literatures [9] . An 
advantage starting from the LW functional is in the con- 
servation laws; they apparently hold if our approximated 
£l[g] keeps the corresponding symmetries; such approxi- 
mations are called as the Baym-Kadanoff approximation 
[10]. Our approximation for .E xc [(/] is 



E xc [g] 



El? A \n\ 
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where the SF functional fl SF contains just the 
temperature-dependent part of the SF energy written as 
ftSF = n SF - n SF0 , where fi SF0 is fi SF at T = so as 
to remove the temperature- independent part. We use 
E^ A [n] given by [11]. This Q SF should have strong 
temperature-dependence. The stationary-point condi- 
tion Sn[g]/Sg = gives g CT (l, 2); it is written as 



Vlda W = V£t M + ^cou (r) + ^lda_*c W , 
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g cr {r l -v 2 ,uJ p ) = -8(v l -Y2), (5) 



where a = ±1 for | and j. <5S^ SF /<5<7 is the contribu- 
tion to the self-energy from SF. w p denotes the fermion 
Matsubara frequencies. 

It is difficult to include the full non-locality and in- 
dependence of 5fl SF /Sg. So we restrict the variational 
space as follows and try to find the optimum solution 
within it. Vlda gives exact solutions of f2[g] if we omit 
Q SF . Therefore we take a variational space of <7 CT (1,2) 
so that it is generated by Vlda plus the onsite magnetic 
fields B R just acting on only d electrons in each atomic- 
sphere (AS) R. Thus g a (l, 2) within the variational space 
is generated by solving a set of SC Eqs. for given B R . 
The set consists of Eq.(5) and 

Kff(ri-r 2 ) = V L V(ri)^(ri - r 2 ) +^^^4 (6) 



Here P R = 6(n - r 2 ) E ro ^(fii, 0i)*2* m (fe, 02) is the 
projection operator to the onsite channel. Note that 
the set {-Br} for all i? is only the variational parame- 
ters, which determine g a (l, 2). Therefore, our problem is 
reduced to searching stationary point of Q[g] as a func- 
tional {Br}. (As in our static approximation E xc [g] ss 
E^ A [n], the stationary point means minimum.) This 
means that we have to solve = 5fl/Sg x Sg/5B R = 

(-V eS + Vlda + 5Q SF /5g^ x Sg/8B R . This determines 
Br as the solution of an inversion equation 



§n sF 5g _ dMR , SF 
Sg sb r 2-? dB R R " 
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where M fl = ^ B f R d 3 r(n^ - nf). B R = B 
if we add external magnetic field B R Kt as a probe. M R 
is the spin magnetic moment of d electrons in R. This 
kind of inversion equation like Eq.(7) is also used in 
the optimized-effective-potential method [12]. Equations 
(5), (6), and (7) give a set of SC equations. With the 
restriction of the variational space, all the quantities in 
RHS in Eq.(l) are treated as the functional of B R in 
Eq.(6), or corresponding magnetic field M R . 

As for f2 SF [g], we give it in a simple model in the on- 
site approximation; it consists of contributions from each 
sites as fi SF [.g] = J2 R ^WISr]- Here g R denotes the d- 
channel onsite part of g a (l,2) as 



g R (r, r', w p ) = .CkKW 7- , ^p) 

m 
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Here we now take just the diagonal part with respect 
to the magnetic quantum number m. &(r,iu> p ) denotes 
solutions of the radial Schrodinger equation. Q R F [g R ] is 
given as 



+|E^ RPA (^x « T K)), 
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(s) = log(l - s) + s 
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v q denotes the bosonic Matsubara frequencies. The 
bare longitudinal and transversal onsite spin-polarization 
functions, X° R (^q) = (S z rS zR ) and x R {v q ) = 
t;(S xR S xR + S yR S yR )o, are given as 

Xfl K) = 7^ WmR&pWmRiVq + Up) 

+9LR( UJ p)g l mR {v q +uj p )}W mR (iuj p ,w q ), (11) 

X° R {v q ) = Yl 2 {dlnR^dLR^n + Wp) 
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where W mR is the square of the radial part of integrals. 
iJ max denotes the radius of AS R. Q R F [g R ] is the SF en- 
ergy given by the onsite-only RPA, though we treat the 
effective interactions I R and I R as the functional of B R . 
They are determined in a SC manner explained later. 
Together with ft R F0 calculated in the same procedure 
with taking the T — > limit in Eq.(9), we can obtain 

^[9] = E R W> wherc = ^R - »| F °- Corre- 

sponding to our onsite approximation on £1 SF [g] , we only 
take diagonal terms on R in the inversion Eq.(7). That 
is, we just take 5g R /SB R , and dM R /dB R , neglecting the 
off site contributions. Then Eq.(7) is simplified as 



oSF _ ( Sg R \ (dM R \ 1 
R \ Sg R 5B R \dB R ) ■ 



(14) 



In practice, X R { v q)i which denotes each term contained 
in x° R and x a R of Eqs. (11, 12), is calculated through these 
formulas symbolically written as 
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T°{v)dv 
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x°M = 

r°(i/)= ldu{f(u)-f{u+w)}X(v+u)Y{u)W mR (u,v). (16) 



Here f(u) in this convolution Eq.(16) denotes the Fermi 
distribution function; X(w),Y(w) denotes the partial 
density of states (DOS) D^ r (lu) = -hng^^/n. The 
factor W mR (Ld, v) introduces a natural cutoff for the con- 
volution of DOS in Eq.(16). 

It is necessary to include contributions due to 
(S™ R S™ R )o which carries SF with other magnetic angu- 
lar momentum m ^ in addition to the contribution 
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of m = (S" l R =0 is S z r given above). If we assume 
spherical symmetry, we can evaluate it as {S™rS™r)o « 
o-\m\ (S zR S zR } Q by counting the number of allowed com- 



binations of mi and mi for m = m\ 
include the contributions by using 



iRPA 
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(17) 



instead of Qq FA (s) in Eq.(lO). Our system is not spher- 
ically symmetric but we utilize this just as a convenient 
method in order to evaluate tt R F [gn] only through quan- 
tities Xr L an d X°r m addition to I R and I R . 

A key point in our method is in the SC determina- 
tion of I R and I R . Here we only treat the paramag- 
netic ground state at finite-temperature. At first, let us 
consider the case without f2 SF [g]; this means that Cl[g] 
in Eq.(l) is reduced to be the crude finite-temperature 
LDA. Then D^^uj) in each site R is given as a functional 
of B R as the solution of an impurity problem; adding 
the perturbation B R xt only at a site R. This problem 
is easily treated in AkaiKKR. Then x5j L (0) is calculated 
through Eqs.(15,16) from D° nR as a function of B c R l . 
On the other hand, we can also calculate the longitu- 
dinal local static spin susceptibility at the site R from 
the numerical 2nd derivative of Cl[g] with respect to B R xt 



Xr(0) = (S zR S zR ) 



1 BMr 
4/4 dB c * 



Note that these 



X% (0) and x^(0) are calculated as functions of B^ xt . At 



1 Mr 



and 



the same time, we also calculate Yr(0) - 

ARy ' 4/4 B c R Kt 

X?? T (fj); we evaluate X/j(0) as if our system is spherical 

symmetric for simplicity. Then we can calculate if so 

that (x^(O)r 1 = {xfWr 1 - if- With this if, we 
can calculate the longitudinal part of fifj F . The same 
method is used also for if so as to be {xh(0)} 1 = 

{x°r{0)} ~ - Ir T - As for fi| F0 , we can calculate it in 
the same manner above from the same D^^uj), Xr{0), 
and x5(0), but taking the T — > limit in Eq.(9) through 
the Matsubara frequencies. As a result, we can obtain 
OfF(M fl ) as a function of Mr. This nf(¥ fl ) should 
give an effect to make the system paramagnetic above 
T c . We should have to include the effect from the be- 
ginning, that is, we rather have to start from Cl[g] in- 
cluding TW[g] = J2r^W( m r)- This means that we 
can obtain a new function Q R (Mr) if we start from 

£l[g] with a trial function D, r (Mr). We impose a self- 
consistency condition, namely the agreement between 
these functions. In practice, we assume a simple form 
fi| F = f^ F + ia SF (Mij/^ B ) 2 - Then our self-consistency 
is just for the parameter a SF , as the trial af F and the cal- 
culated aQ F t are the same. This assumption is necessary 
so as to make the SC equation determining Q r f (Mr) be 
closed. We can take the self-consistency used here as a 



simplified version of the general self-consistency scheme 
in the construction of the LW functional; the 2nd deriva- 
tive of il[g] with respect to g a (l,2) can give the fluctu- 
ation (or 2-body propagators); on the other hand, the 
fluctuation can be used to construct r2 X c[ff] through the 
coupling-constant-integral formula [12]. In the STLS the- 
ory [13] for homogeneous electron gas and the SCR the- 
ory [5], other simplified versions of the general scheme 
were used. 

FIG. 1. Determination of a SF : This is an example for Fe. 
We determine a SF so that af F — a^ t .See text. 
0.03r 
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In Fig.l, we give an example of plots of af n F vs. a^ t , 
determining a SF . In our formalism, the interaction if 
has a part that is linearly-dependent on af n F . The most 
divergent part in Eq.(9) occurs at vq = in log(l — 

ifxfiO)) = log^WrV^O)); log(x5i(0)) = 

- log (x R (0) + 4of n F ) , where X r(0) denotes X r(0) at 

ofF = 0. It diverges at af n F -> -\ {x^( )} • This 

occurs at ~ af F = 0.0139 Ry denoted with the arrow 
in Fig.l. In the case of Fe shown in Fig.l and fcc-Co, 

— |xfl(0)| is positive even at T = because LDA 
gives a local static magnetic moment [16]. Therefore our 
treatment gives a SF -> -\ |x«(0)} > at T -> 0, 

thus Xfl(0) — > oo. Thus fi SF keeps Xr(0) > and for- 
bids the existence of the static local magnetic moment at 
T = 0. This is physically reasonable as a ground state 
character, though we still have an unphysical divergence 
at T — » 0. In these cases, the results does not coincide 
with the LDA limit. On the other hand, as for Ni and 

Pd, which show negative — |x#(0) j , it coincides with 

the LDA limit with a SF -» at T -» 0. 
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FIG. 2. Top Panel: Bulk Magnetic susceptibility 
The middle panel: The local onsite susceptibilities 
{Xfl(0)} _1 > {x I r j (0)}~ 1 , and the effective onsite interaction 
ft = {Xh(O)}" 1 - {x'r (0)}" 1 . The loser panel: a SF in 



1 ft 

qSF _ Q SF 



or (Mk//Ub) ■ Solids lines for experimental 
lattice constants, broken lines for 3% smaller lattice constants. 
See text. 
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The partial DOS D^ r (uj) are calculated with ~ 8000 
k point in the irreducible Brillouin zone (BZ), where 
each poles of eigenvalues are effectively smeared with 
Lorentzian functions oc l/(^ 2 + A 2 ). We needed to per- 
form calculations with A L = 0.009, 0.012, 0.016, 0.022Ry 
cases and to take an extrapolation Al — ► 0. After we 
obtain the self-consistently determined a SF , we can cal- 
culate the bulk susceptibility % f° r the whole systems. 
In Fig. 2, we summarise our results. Used lattice con- 
stants a = 5.47, 6.70, and 6.80 a.u. for bcc Fe, fee Ni, 
and fcc-Co (experimental values at T c ); a =7.35 a.u. for 
fee Pd (experimental value at T = 0). They are shown 
with solid lines. In addition, we show the cases for 3% 
smaller lattice constants by the broken lines; the differ- 
ences from solid lines indicate errors due to LDA be- 
cause the ordinary LDA calculations give rise to errors 
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of this level in the lattice constants. The bottom panels 
show a SF . The top panels are for the inverse of the bulk 
magnetic susceptibilities x -1 calculated from these a SF , 
where we also show experimental results taken from Ref. 
[17]. The middle panels are for {x?j L (0)} _1 ,{xk(0)} _1 
and/={x° fl L (0)}- 1 -{x^(0)}" 1 . 

As for Pd, our LDA calculations without fl SF in ASA 
at T = give ferromagnetism for a— 7.35 a.u. and para- 
magnetism for a = 7.13 a.u. In Ref. [14], the LDA calcu- 
lation with a — 7.35 gives paramagnetic phase. This dis- 
crepancy might arise from the difference in the formalism 
and details of implementation [15]. However, the differ- 
ence seems not important for the overall discussions here. 
The temperature dependence of x _1 , especially slope at a 
high temperature region, is rather well described in spite 
of no tunable parameters. It fails to describe x _1 around 
T = where we expect only SF around k ~ dominates 
the behaviour. It is reasonable because our formalism 
could work well only in the case that SF in each site be- 
haves independently; in other words, all the modes of SF 
are excited rather homogeneously in BZ. As T becomes 
higher, a SF increases smoothly, which suppresses SF at 
k = (x _1 gets larger). On the other hand, the local SF 
Xfl'(O) does not change much. The enhancement factor 
on local SF Xr(0)/x'r j (0) due to the interaction 1^, is not 
so large <~ 1.5 at T = and gradually decreases for larger 
T. In Ni, our method can also describe experimental T c 
and the slope of x 1 reasonably a SF decreases above 
T ~ 1000K. x" 1 is n °t so linear as experiments and 
contains negative-quadratic dependence as x _1 ~T 2 
reflecting the behaviour of a SF . In both cases of Pd and 
Ni, it seems that we have to add positive-quadratic de- 
pendence x _1 oc T 2 , so as to have better agreement with 
experiments. We might expect such an effect if we in- 
clude k-dependent SF, because then we will have larger 
contributions even for lower temperatures. 

As for Fe and fcc-Co, calculated x _1 is rather linear. 
Slopes of x are too large compared with experiments, 
though calculated T c 's are still reasonable. The factors 
Xk(0)/x?j j (0) are rather large ~ 5 in these cases. The re- 
lation {xk(°)}~ at T should be hold m agree- 
ment with the previous discussion. This condition keeps 
the temperature dependence of (x#(0)} too large, cor- 
responding to too large a slope of x ■ As is same the 
as cases in Ni and Pd, we expect better agreements by 
taking account of k-dependent SF which will remove the 
above condition from our formalism. 

We have shown results of paramagnetic susceptibilities 
for some metals. They give reasonable agreements with 
experiments, indicating the possibilities of treatments 
along the present line including SF self-consistently. In 
particular, the results seems to indicate necessity of in- 
cluding the k-dependence of SF for the further improve- 
ment. 

We thank Dr.Aryasetiawan for valuable discussions 



and comments on the manuscript. 
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